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We present a clear and rigorous derivation of the Ewald-like method for calculation of the 
electrostatic energy of the systems infinitely periodic in two-dimensions and of finite size in the 
, third dimension (slabs). We have generalized this method originally developed by Rhee et al. [Phys. 

Rev. B 40, 36(1989)] to account for charge-dipole and dipole-dipole interactions and therefore made 
it suitable for treatment of polarizable systems. This method has the advantage over exact methods 
of being significantly faster and therefore appropriate for large-scale molecular dynamics simulations. 
It however involves a Taylor expansion which has to be demonstrated to be of sufficient order. The 
method was extensively benchmarked against the exact methods by Leckner and Parry. We found 
it necessary to increase the order of the multipole expansion from 4 (as in original work by Rhee et 
al.) to 6. In this case the method is adequate for aspect ratios (thickness/shortest side length of 
the unit cell) < 0.5. Molecular dynamics simulations using the transferable/polarizable model by 
Rustad et al. were applied to study the surface relaxation of the nonhydroxylated, hydroxylated, 
and solvated surfaces of a-Fe2C>3 (hematite). We find that our nonhydroxylated structures and 
energies are in good agreement with previous LDA calculations on a- alumina by Manassidis et al. 
[Surf. Sci. Lett. 285, L517, 1993]. Using the results of molecular dynamics simulations of solvated 
interfaces, we define end-member hydroxylated-hydrated states for the surfaces which are used in 
energy minimization calculations. We find that hydration has a small effect on the surface structure, 
but that hydroxylation has a significant effect. Our calculations, both for gas-phase and solution- 
phase adsorption, predict a greater amount of hydroxylation for the a-Fe2C>3 (012) surface than for 
the (001) surface. Our simulations also indicate the presence of four-fold coordinated iron ions on 
the (001) surface. 
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I. INTRODUCTION 
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In a wide variety of scientific and engineering problems, simulation techniques connecting models at Angstrom 
scales to 'microstructural' scales 10 4 -10 5 times larger are required. Classical simulation of atomic and molecular 
systems works between these two extremes and is useful for moving atomic scale predictions into the microstructural 
or thermodynamic regime . 

Here, we consider simulations using models with analytical interaction potentials involving charges and dipoles in 
interfacial environments. The long-range coulomb interactions in such systems must typically be taken into account 
over distances much larger than the sizes of samples which can be directly simulated. Historically this problem has been 
i-^j ■ approached by infinitely and periodically repeating the simulation sample in all three dimensions. The long-ranged 
contributions to the energies and forces can be determined by a variety of methods. They include transformation to 
rapidly convergent sums in real space [|J , the fast multipole method |s) and various Ewald-type methods @^-§. The 
latter involve summations over the reciprocal lattice of the 2-D lattice. The relative merits of different methods are 
reviewed by Tildesley JlOj ] . The original 2-D Ewald method due to Parry Q is commonly used in quantum mechanical 
calculations where the Ewald summation comprises a small fraction of the computation time. For simulations with 
parameterized analytical potentials, using Parry's method makes the Ewald summation procedure rate-limiting and 
it becomes worthwhile to explore alternatives to Parry's method. 

In the present paper we extend the 2-D periodic Ewald method of reference Q to include polarizable atomic 
centers. Effective pairwise potential energy functions parameterized for bulk systems do not work well at surfaces 
because the surface intermolecular environment is entirely different from that of the bulk. The first step in enhancing 
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the transferability of empirical potentials is to consider explicitly the effects of electronic polarization [[Ll). Therefore, 
the extension to polarizable systems is quite important in a method designed to be used at interfaces. In addition, 
we have fixed some omitted fourth order terms in reference || and we have carried out the Taylor expansion to sixth 
order in the aspect ratio, significantly improving the accuracy of the method. We have also clarified the extension to 
charged systems. 

The generalized method is tested against previously-published benchmarks, and then applied to investigate the 
structure and hydroxylation of a-Fe203 (hematite) surfaces. Hematite is chosen because it (as well as its aluminum 
analog, corundum) is a good example of a material which is well-studied using high-vacuum techniques and state of 
the art theoretical methods, and is also important in natural weathering and adsorption processes of interest in soil 
chemistry fl2| . 



II. METHODS 



Our polarizable system is most conveniently treated by generalizing the expressions for charge-charge interactions 
from Ref. Q to include charge-dipole and dipole-dipole energies and forces. 

We start by deriving the Ewald-type expression for the energy of a system of point charges periodic in two dimensions 
and finite along the third dimension. If the monopole-monopole potential energy expression is given by < i ) mm (r, s, z), 
then one can work out the expressions for the charge-dipole, and dipole-dipole interactions using: 

$ md (r, s, z) = -qjili ■ V$ mm (r, s, z) (1) 

for the charge-dipole interaction, and 

$dd{r, s, z) = -(pj ■ V)(|Ui • V)$ mm (r, s, z) (2) 

for the dipole-dipole interaction, as well as gradients involving these terms arising in the expressions for the forces 
& 

The derivation of the Ewald summation for the case of two-dimensional slabs of finite thickness is given in previous 
publications || and ||. We present another method of derivation that we believe is easier to follow and rigorous. It 
is essentially a special case of the derivation presented in H . 

The electrostatic energy of the system of point charges with 2-D periodicity is given by 

*™^EE'^4<X>> J i/f (3) 

Here fy is the separation between particles i and j, Sij is the projection of this vector on the plane, v is the translation 
vector in the plane. N is the total number of particles, and A is the surface of the 2-D unit cell. The last term is 
the energy of the interaction of the slab (if non-neutral) with the uniform neutralizing background. The prime on the 
summation sign means that the term i = j is skipped when v = 0. A vector quantity without an arrow denotes the 
magnitude of the vector. 

As in the 3-D Ewald method, we split the infinite sum (||) into short-ranged (Vs) and long-ranged (Vl) parts. Vl 
will be replaced by its Fourier series expansion. To perform this separation, we write the identity 

1 = 1 _ ho(s,j) + fop(s,7) ^ 

r r s s ' 

where ho(s,j) is the convergence function with a parameter 7 (to be defined below). One then obtains: 

*mm = V S + V L , (5) 

v s = \j2J2' qiq i i 1 l r v v - h o{ s ijP,l)/ s ij») ( 6 ) 
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The cutoff function h(s)/s is denned by specifying its Fourier image. Consider the function w(s) = ho( s v)/ s v, 
where s v = s + v and s v — |s„|. It is periodic with respect to the 2-D translation vector v. Therefore it can be 
expanded in Fourier series in 2-D reciprocal space vector G, being the Fourier coefficients: 

w(s) ee 2J h Q {s v )/s y = 2J Cq exp(iG ■ s) (8) 

The sum is over all reciprocal lattice vectors, including G = 0. Following the same notation as in Ref. ||, we define 
the function /(G) from 

where A is the area of the 2-D unit cell. The Fourier coefficients Cg of the function J2ffho( s v)/sv are 

C g = jJ a J2 ^ ex P (-zG • 3)dS (10) 
Changing the order of summation and integration and using the additivity of integration, one can write 

C g = t/ expHG ■ s)ds (11) 

J all s— space " 

Comparing the expression above to the definition of the 2-D Fourier transformation of the function ho(s)/s, 



/ exp(— zG • sjas, 

J all s — space " 



2tt 

one obtains 

f(G) 1 f h (s) 



expHG • s)ds. (12) 

" Jail s— space s 

This means that /(G)/G is the 2-D Fourier image of ho(s)/s. Following Ref. H], we choose 

/(G) = exp[-( 7 G) 3 ]. (13) 
The direct space cutoff function ho(s)/s is obtained by transforming back from G-space to s-space 

ho{s) - 1 f eMiG-^ M{ ^ Gf] dG (14) 



2 71 " Jail G-spaco G 

This integration is conveniently done using polar coordinates in the G-plane: 

exp[(- 7 G) 3 



_2vr 

Here we used the well-known integral 



exp(iGs cos 6)d6 



Jo(Gs)exp(- 7 3 G j )dG. (15) 



exp(ix cos 9)d8 = Jo{x), (16) 



where Jq(x) is the Bessel function. Vl will be calculated as the sum of Fourier series. To do so, rewrite it using the 
periodic function w(s) defined above. This is done by adding and subtracting i — j term. 



n - iEE^(»«) - 5 (a^) Erf - n (f>) L 



i j l — l \ l — l 



- (") 

all s— space 
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The limit in the expression above exists, and the direct calculation yields Q 



s->o s 37 \ 3 
Substituting (|l8|) and (111) in (|l7|), one obtains 



limMf) = J_r(T). (18) 



- (19) 

all s— space 



a „ fj i=l / \ / i=1 / J a. 



all G 

This expression contains two infinite terms, the G = term in the Fourier sum and the last term. The sum of these 

r h (s)- 

-/all s — space s 



terms, I = f ,, fe oW — I. fig can k e written as 

' Jail s — space s ' 



I = lim 

G^O 



1 /" ho(s) , .a . If , .a .ds 

— I exp(— iG ■ s)ds — - / exp(— iG ■ s) — 

^ J all s— space ^ Jail s — space ° . 



(20) 



As one can see from the Eqs. (||) and ( ||), the first integral in ( pp|) is ^J-^j • The second integral is calculated by 
transforming to polar coordinates and using the formula 

Jo(bx)dx — 1/6. 

Its value is jfe. Combining these expressions, one obtains: 

I = lim 2 ^ m ~ l) = hm ^M-^f)-D = Q } 
g^o A G g^o AG y ; 

Thus for our choice of f(G) — exp[(— 7G) 3 ], as well as other functions /(G) satisfying the above condition, the term 
due to the background charge cancels exactly the G = term of the reciprocal space sum. Therefore the long-ranged 
part of the interactions can be written as 

1 = jE 22^1 1 . exp(*3 - *)|» - £r (i) £ „? (22, 

G#0 i=1 ' 



i=l 



Here we used the identity 



N 



giQj exp(iG ■ Sjj) = ^2 1i exp(iG ■ s^l 2 (23) 

i,j i=l 

In practical calculations, the summations in real space are truncated at \v\ = t? max . The cutoff function h(s)/s 
does not dissipate at large distances fast enough like it does for the 3-D Ewald method. Instead it has the asymptotic 
behavior h(s)/s m 1/s at large distances. It is also convenient to separate the i = j term in real space. We therefore 
write 

* mm (r,a,z) = Vt + V L + V m (24) 

W-5 £ «/El|^-MW/^)4E^E(;-^) (») 

£7=0 7 i i/>0 v 7 

f— i i : f iw^rw^) (26) 



1 r(|) = 2.67893853471 
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The expression for Vi is given by Eq. ([22|) . 

The series for V m (Eq. |2g) converges slowly, and direct calculation is extremely inefficient. Therefore it is advanta- 
geous to be able to precompute certain terms. Defining the function v' m 



Ar + i>\ \s + u\ 



one can write 



V™ = \ E (28) 
The key to the method of reference H is the use of a multipole expansion to compute V m . The Taylor expansion of 

(l+x) 1 / 2 ^ 

1 x 3 9 15 , 105 4 

= 1 - - + -x 2 - —x 3 + —x A + ... (29) 



(l + x) 1 / 2 2 8 48 384 

This, in conjunction with the elementary definitions of r and s 

r = [\s + v\ 2 + z 2 ) 1 ' 2 



1/2 

"( ~2 + '- + 2l: Y- + 1 ) ( 3 °) 



s = (\s+pn 



2U/2 



2 - - \ 1/2 



I/ 2 V 



and noting that any odd powers of s ■ v will vanish in the sum on v, yields the following expression for the multipole 
expansion in powers of (s,z)/u: 



-, OO-. ~ oc -. ^ oo-. 1r . oo / -* -*s 2 

-Z > — T + -Z S > — + -2 > — k —2 > = 

2 ^ i/ 3 4 ^ j/ 5 8 ^ v b 4 ^ i/ 7 



lg ^z +z j ^ „ 7 48 z 2^ „7 + g 2^ ^9 

L/ — Z^ max zy— z/ max y — y max 



105„ 4 ^ (*.*) 2 
16 ' 



£ ^ (32) 



where we retain only the terms s n z m for which n+m < 6. 

Application of equations [l] and || to obtain $ mc ; and <&dd is straightforward for V s (equation ^) and Vl (equation H|) , 
following closely the 3-D periodic methods jl3) . Reference |Hj gives a generalization of the method to polarizable 
systems. 

To implement these equations we split the interval of interest in 3 parts and approximated ho(s)/s using the ratio 
of two polynomials. The 3 derivatives of ho(s)/s were approximated separately, also using rational approximation. 
The optimal rational approximation can be easily constructed using the code given in Numerical Recipes [fL5| . 

We report here the derivatives of V m (equation |32| ) in the Appendix, as these terms, while straightforward, are 
somewhat tedious to compute. We also note, as pointed out in reference Q, that the expressions given in || do not 
contain all the terms to the fourth order. The expressions in the Appendix include multipole contributions up to the 
sixth order. Also, the numerical value of Yla.il ^V^ 3 f° r ^ ne sc l uare lattice with parameter L as given in reference 
11 has only 3 correct digits and should be replaced with 9.0333388/L 3 . This result was calculated with quadruple 
precision. 
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III. BENCHMARKS AND CONVERGENCE TEST 



As one can see from the previous section, the Ewald-type method for 2-D slabs has 4 computational parameters: 
the convergence parameter 7, cut-off parameters for summations in real and reciprocal space and the order of Taylor 
expansion used to calculate V m . In this study we used 7 = 0.1L s i lort , where L s h or t is the shortest side of the 
rectangular unit cell. The summations in real space included the first shell of neighboring cells, the total number of 
cells included being 9. We increased the order of expansion from 4 (as was originally used by Rhee et al § and also 
used in Ref. &) to 6, and this improved the accuracy of the calculations. 

A valid test of convergence with respect to all computational parameters is a benchmark against another method. 
A test example for the monopole-monopole interaction was provided in Ref. ]10| . As discussed in that review, the 
method of Leckner Q has only one computational parameter, the cutoff of the summation in real space. It has 
exponential convergence and can be made as exact as required by choosing this cutoff appropriately. Therefore it is 
ideal as a benchmark Jig] , The comparison of the forces and potential energies calculated in this study with those 
from Ref. jl(^] and calculated using the implementation of Parry's method || is given in Tables [j] through III . As one 
can see from Tables | and ||, the two benchmark methods (Leckner's and Parry's) yield exactly the same results. One 
can therefore be confident in the benchmark itself. The results of our method are in excellent agreement with the 
benchmark for small interparticle separations in z direction. Once the ratio of the maximum separation in z direction 
to the magnitude of the translation vector (aspect ratio of the slab) increases, the agreement degrades slightly, the 
maximum error being 2.5 % for the aspect ratio 0.5. The Taylor expansion used here to efficiently evaluate V m is 
in powers of {zij/v) 2 . Therefore the quality of representation of V m degrades as the expansion parameter increases. 
This explains the results of the benchmark. The error of 2.5 % in F z can be considered acceptable. Real molecular 
dynamics simulations are done with ^1000 atoms, not with 2 atoms. As noted by Tildesley |10|, the convergence of 
2-D Ewald-type methods improves when the cell is filled with charges, so that ions on opposite faces of the slab are 
shielded from each other. 

We performed a test on the slab of 216 formula units of Fe2C>3 (1080 atoms) that corresponds to the neutral 
termination of the (012) face. The electrostatic part of the total energy and forces in this system were calculated 
using Parry's method and the method described in this study. The calculations yielded the total energy of -2584.148 
e 2 /A and -2583.626 e 2 /A respectively. The maximum deviation in the magnitude of any force on any atom was 1.17 
% of the exact value. 

To the best of our knowledge, the benchmarks for charge-dipole and dipole-dipole interactions in 2-D periodic 
boundary conditions have not been published. We therefore present them in Tables |v| and [v[ respectively. Once 
dipoles are present, the energies calculated using the Parry's method Q and our method differ by an additive constant. 
This stems from the fact that the surface terms arc treated differently. The agreement of forces is good. Since evalu- 
ation of charge-dipole and dipole-dipole interactions involves successive differentiation of V m , the order of expansion 
decreases. Therefore it is not surprising that the quality of calculated charge-dipole and dipole-dipole forces degrades 
faster with the aspect ratio of the slab than does the force for the case of charge-charge interactions. Overall, the 
agreement is still totally satisfactory. Since the agreement with an independent method was good, we conclude that 
the choice of the computational parameters is adequate. 



IV. APPLICATION TO HEMATITE SURFACES: BACKGROUND AND PREVIOUS WORK 

The morphology of hematite (a-Fe20s) crystals is complex p^-pl[, both in natural and synthetic samples. The wide 
variety of crystal faces present has tended to discourage atomistic interpretation of binding at hematite surfaces, with 
some notable exceptions [ ^6|Jl7| ]. This contrasts with goethite (a-FeOOH), where there have been several experiments 
and much discussion concerning the role of contiguous surface (monocoordinated) hydroxo groups in bidentate binding 
of oxyanions p8| , p2^ |. 

The properties of both hematite (a-Fe203) and corundum (a-A^O^) surfaces in vacuo have been studied experimen- 
tally p^Jl9| , p4] -p6[ and by computer simulations p7j-p3|. In this section we consider both ideal vacuum-terminated 
and hydroxylated (001) and (012) surfaces of hematite. These are by far the best studied, and among the most 
prevalent surfaces for both hematite and corundum |24||34| . 

As is the case for any oxide surface, neither the extent of hydroxylation of hematite, nor its hydroxylated structure is 
well-understood p5| . Surface hydroxylation has been demonstrated to occur on hematite by vibrational spectroscopy 
jl8],[36|], but the spectra are very complicated and not easily interpreted. Here, we offer an application of our 2-D 
periodic Ewald summation method for polarizable systems to investigate the surface structures using our previously- 
developed methods for treating Fe-OH and Fe-OH 2 interactions P,p7|. 
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V. APPLICATION TO HEMATITE SURFACES: SIMULATION PROCEDURE 



The atomic interactions were described using parameterized potential functions developed previously |3q ]. Briefly, 
this model includes formal charges on all ionic species, polarizable oxygen anions with a cut-off function to reduce 
induced polarization at short interionic separations, and various two- and three-body short-ranged terms. The func- 
tional form of all interactions and the list of parameters is given in Ref. |Q . The lattice parameters of hematite at 
zero pressure were determined in Ref. |j38j by minimizing the total energy of the system with respect to the magnitudes 
of lattice vectors and internal structural parameters. They are listed in Table 3 of Ref. |38[] along with experimental 
values. It is worth pointing out that the 0-0 and H-0 interaction parameters used here were exactly the same as in 
the model for liquid water and gas phase clusters of water molecules. Considering this stringent transferability require- 
ment the agreement of the calculated lattice constants with experimental values can be considered good. However, 
the Fe-0 bonds predicted by the model are systematically about 3%-4% too long. 

The simulations were done in two different geometries: in "sandwich" geometry with 3-D periodic boundary con- 
ditions (PBC) and in slab geometry. In simulations with 3-D PBC the length of the translational vector along z 
direction was 3 nm for all faces. Therefore the slabs were separated from each other by empty space in z-direction 
(~1.5 nm thick). 

The atomic positions were relaxed by the following procedure: The sample was heated to 400 K and cooled to K 
in 3-D molecular dynamics (MD) simulations 2 ps long. Then the energy of the system was minimized with respect 
to atom positions using the individual damping method by Bennett fl39| |. In this method, the equations of motion are 
solved as in standard MD but the velocity Vi of an atom is set to zero if Fi ■ Vi < 0, Fi being the force acting on the 
atom. This procedure is carried out until convergence of the potential energy is reached. The 3-D minimal structure 
was used for the initial configuration for 2-D energy minimization (with no periodicity along z-direction). 

Our surface energies are calculated using 

l = \{E slab -E bulk )/A (33) 

where E s i ab is the potential energy of the slab system, E{, u ik is the potential energy of an equivalent number of Fe203 
molecules in the bulk, and A is the surface area. 



VI. RESULTS 



A. Nonhydroxylated Surfaces 

The main utility of our potential set is our ability to treat surface hydroxylation. However, because there is little 
basis for comparison with other theoretical or experimental results on hydroxylated surfaces, we are compelled to test 
the performance of our potential set on non-hydroxylated surfaces. Even for these, there is little information available 
on hematite specifically. However, as mentioned above, there has been a substantial amount of work on corundum 



surfaces to which we can qualitatively compare our results [f29|-pl 33 32y . We limit the investigation to the structure 
and energetics of the (001) and (012) surfaces, as these are the surfaces receive the most attention both theoretically 
and experimentally. 



1. (012) Surface 



The (012) surface has been previously studied in references [ p4|40]j25| , |30| . The (012) surface is often indexed using 
the symmetrically equivalent (1102) surface in hexagonal notation; also note that in reference p3] this surface is 



is shown in Figure 0(a) • This surface can 
. The slab in Figure 0(a) is comprised of 



indexed as (1012). A cross section of the (012) slab, as viewed down [100 
be built up from neutral O-Fe-O-Fe-0 units as described in reference |}0| 
three such stacking units. The unit cell predicted by the model is 0.556 nm x 0.518 nm. For the calculations here, 
we use a supercell consisting of 6 x 6 unit cells; the supercell is periodic in X and Y directions. 

The atomic coordinates for the unit cell of the slab are given in Table [v]. In Figure | we show the two side views 



and the map view of the relaxed (012) top stacking unit (including the top O layer from the middle stacking unit). 
The noteworthy zig-zag oxygen ridges characterizing the (012) surface are clearly seen in Figure 0(b) . The "slopes" of 
the oxygen ridges are occupied by five-fold coordinated iron ions. These are the darkly shaded iron ions in Figure 0. 
To test convergence of slab energies and structures with respect to slab thickness, we considered systems having one, 
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two, three, and four O-Fe-O-Fe-0 layers. From Tables VII and VIII, one concludes that 3 layers is sufficient for 
convergence for both energy and structure. The thickness of the nonhydroxylated 3 layer slab, shown in Figure |l| is 
0.964 nm, and, for the size of the computational supercell considered here, the system consists of 216 Fe20a_molecules 



(1,080 total atoms). The relaxation of each layer in the direction normal to the surface is shown in Table VII. Upon 
relaxation, the top three layers are brought closer together, compensated by separation of the bottom O-Fe layer. 
The Fe-0 distance across this layer increases from 216 pm to 228 pm. 

There is, to our knowledge, only one theoretical study of a alumina (012) |50|. That work employed embedded- 
cluster LDA me thods to calculate energies of unrelaxed surfaces. Our size-converged unrelaxed surface energy of 2.8 
J/m 2 (see Table VIII) compares well with the value of 2.95 J/m 2 obtained there. We find that allowing for relaxation 
decreases the surface energy to 2.0 J/m 2 . 



2. (001 ) Surface 

The (001) surfaces of both hematite and corundum have been studied extensively J24|,[ll|]. We consider here the 
standard neutral, nonpolar, nonhydroxylated termination based on the stacking of Fe-O-Fe building blocks, as de- 



scribed in |30|. A cross section of the four Fe-O-Fe layer (001) slab, viewed down [100], is shown in Figure |2j(a). The 
primitive hexagonal unit cell has a=0.518 nm; here we work with a nonprimitive orthorhombic cell with dimension 
a=0.898 nm and b=0.518 nm. The computational supercell is 4x6 orthorhombic unit cells. An orthorhombic cell was 
chosen for computational efficiency. 



The atomic coordinates for the unit cell of the slab are given in Table IX. In Figure |2|, we show the two side 
views and the map view of the relaxed (001) surface. We include the top two Fe-O-Fe stacking units, but omit 
the bottom Fe layer for clarity. In other words, we show the top octahedral layer plus the overlying layer of iron 
ions. The model surface structure is characterized by a top layer three-fold coordinated iron ions (darkly shaded in 
Figure]^) in hexagonal arrangement. Configurations with two, three, four, five, and six such blocks were used to test 
for convergence of surface energy and surface structure as a function of the thickness of the slab. Table [x| shows the 
percent change in distances between layers parallel to the surface with respect to bulk distances. Relaxation brings 
both the top Fe and O layers and the two interior Fe layers closer together. 



As shown in Table XI , our calculations are in good agreement with previous calculations on isostructural a alumina, 
|^,|l]j33) except that in the tight-binding calculations of reference |33| , the distance between the two intra-octahedral 
layers of Al ions increased by 20%; in contrast to the 40-50 % decreas e pr edicted by other empirical models and LDA 



51[ . The relaxed and unrelaxed surface energies are reported in Table XII . These compare favorably with the energies 
for a alumina reported in reference J3lj , as well as the unrelaxed surface energy reported in pC| (see Table XIII ) . The 
structure, energy, and relaxation energy of the (001) slab appear to be well-converged by the fourth, (or even third) 
Fe-O-Fe layer. The thickness of the four-layer slab, shown in Figure ||(a) is 0.773 nm. For the X and Y dimensions of 
the computational supercell, the four-layer slab consists of 224 molecules (1120 total atoms). The convergence with 
thickness is again similar to what was found for a alumina in reference [ pl| . 

Experimental fracture surface energies, as discussed in reference p0[ , range between 24 J/m 2 and 6 J/m 2 for (012) 
to around 40 J/m 2 for (001). The reason for the uniformly low values estimated by all theoretical methods cannot be 
speculated on at this time, but is a subject for future work. 



3. 2-D Periodic vs. 3-D Periodic Calculations 



As one can see from Tables VIII and XII, the minimal energy for the case of 2-D PBC was always slightly higher 
than that for 3-D PBC. This difference can however be considered small. We have therefore demonstrated convergence 
with respect to the thickness of the empty space for 3-D calculations. The difference in atomic positions between 
3-D and 2-D minimal energy structures was no more than 0.002 nm for all surfaces considered. The fact the 2-D 
calculations agree well with the 3-D (which don't have any Taylor expansion) in case of sufficient separations between 
slabs in z direction suggests that for the actual system of interest the convergence with respect to computational 
parameters was reached. This is of further interest because 3-D PBC need to be implemented in electronic structure 
methods incorporating plane waves |3l| . Our calculations support the inherent assumption that repeating the slab in 
the z direction with sufficient empty space between the layers has only a small effect on surface structures and energies 
for the systems considered here. This conclusion is of course system-dependent and would need to be demonstrated 
in any given case. 
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B. Hydroxylated Surfaces 



In the previous section, it was demonstrated that our interaction potential was in good agreement with other 
theoretical studies for predictions on nonhydroxylated surfaces. Extension to hydroxylated systems is quite difficult 
for both the Hartee-Fock and LDA methods. It has been demonstrated that treatment of hydrogen bonds requires 
very large basis sets and correlated wave functions for Hartree-Fock-based methods Q and gradient corrections are 
required for density-functional methods 43 1. Our potentials are competitive with these methods. For example, we 
obtain the same ordering of conformers of H30(H20)4 + as [H and the same ordering of the OH(H20)3- as [Ef|. Our 
model will provide a semi-quantitative representation of surface hydrogen bonding. 



1. Generation of the Hydroxylated Surface: Molecular Dynamics Calculations 

Defining the hydroxylated surface structure is in itself a difficult problem PqJIql , The coordination numbers of 
the iron atoms and extent of hydroxylation are unknown. Because our potential model is as general as possible, we 
generate the hydroxylated surface structures through molecular simulation. To begin the simulations of the interfacial 
systems, we take the vacuum terminations of the 216 formula unit (012) system (3 layers of O-Fe-O-Fe-O), and 224 
formula unit (001) system (four layers of O-Fe-O-Fe-O) as starting configurations. The vacuum configurations were 
then placed in contact with a layer of water (about 1.5 nm thick) obtained from a simulation of bulk water at 300 
K. Simulations of both the (012) and (001) interfacial systems were run for 50 picoseconds at 300 K. This time was 
deemed sufficiently long; by about 40 picoseconds, the potential energy had ceased drifting downward and began to 
fluctuate about its mean value. 

The surface structures for the (012) and (001) systems are shown in Figures || and |J, respectively. For the (012) 
surface, as might have been expected, each bare iron ion has indeed adsorbed a water molecule from the solvent. In 
several cases, "chemisorption" has occurred, and some of the FesO sites on the surface have protonated. 

For the (001) surface, however, the assumption of six-fold iron coordination does not appear to be strictly valid; 
most iron atoms are four-fold coordinated. Here there is a balance among the surface relaxation forces (acting to pull 
the iron ion inward toward the surface), the forces of solvation (acting to pull the iron ion out towards the solution), 
and the lateral surface hydrogen bonding forces. Those surface iron ions which are 6-fold coordinated have slightly 
"popped out" of the relaxed top layer of iron ions. Because the relaxation energy for this surface is quite large, this 
evidently requires a significant amount of work. Thus the large surface relaxation energy tends to keep the top iron 
layer in a four-fold coordinated state. 

The issue of the coordination state of the top iron layer is important because x-ray adsorption studies of the 
oxidation of metallic iron have often revealed evidence for the presence of four fold coordinated iron oxides in the 
passivation layer [Q. This has been used to suggest that the passivation layer consists of maghemite (7Fe20s). While 
this seems likely due to the additional evidence from x-ray scattering, it is important to note that the presence of 
four-fold coordination of iron does not appear to be a sufficient condition for the presence of maghemite. 

Of course, for a hematite crystal growing in the presence of a supersaturated solution of iron ions, this artificially 
cleaved and hydroxylated surface will adsorb iron ions and grow. This study obviously says nothing about the 
mechanisms of iron adsorption and crystal growth. For example, it is possible that our artificially-prepared surface 
has an extremely high affinity for iron and would never be realized for a surface in equilibrium with iron in solution. 
In other words it may only exist in nature over, say, microsecond time intervals (where no growth occurs) after the 
solvation of a freshly fractured hematite crystal. At the present time, however, it represents our best, least-biased 
guess for the structure of a hydroxylated hematite (001) surface. 

These molecular dynamics studies are meant to reveal the coarse aspects of the surface structure. Interpretation 
beyond this is unwarranted because the thermodynamic states of the systems during hydroxylation are not well- 
defined. In particular, there is no reliable way of defining a solution pH, which is clearly a controlling variable in the 
extent of surface hydroxylation. However, it is possible to make some fairly precise statements about the predictions 
of the model concerning the energetics of hydroxylation at K with no solvent present at the interface. 



2. Structure and Energetics of Gas-phase Hydroxylation: Molecular Statics Calculations 



(012) Hydroxylated Surface Structure 

For the (012) surface, the first step toward "gas phase" surface hydroxylation is the adsorption of a water molecule 
to the five- fold coordinated iron ion (see Figure @) . This produces a chain of water molecules sitting over the "valleys" 







formed at the middle oxygens in the O-Fe-O-Fe-0 sequence, spaced about 300 pm apart. The bond distance to 
the (formerly) five- fold coordinated iron ions is 212 pm. The coo rdinates of atoms in the unit cel l of the slab for 
hydrated and hydroxylated surface are given in Table XIV and XV , respectively. As shown in Table XVI , adsorption 
of water molecules has only a minor effect on the rest of the surface. Hydroxylation of the surface takes place by 
dissociation of the adsorbed water and protonation of the zig-zag ridge of oxygens (Figure ^J) . The rows of adsorbed 
hydroxide ions over the oxygen "valleys" move closer together (280 pm). Hydroxylation has a significant effect on 
the surface structure, moving the top layer of iron ions much closer to the top of the oxygen "ridges" . The bond 
between the top-layer iron ion and the bottom-layer oxygen ion in the O-Fe-O-Fe-0 sequence is elongated from 210 
pm (nonhydrated) to 212 pm (hydrated) to 238 pm (hydroxylated). 

Concerning the speciation of surface hydroxyls, per unit cell, there are two =FeaOH and two =FeiOH sites, with 
a relatively weak 203 pm hydrogen bond forming between the singly and triply coordinated hydroxyls. The 0-0 
distance between =FesOH and =FeiOH is 290 pm. The singly-coordinated sites are nearly colinear along the X 
direction, spaced 280 pm apart. Overall, the arrangement of =FeiOH groups is quite similar to the (110) surface 
of gocthitc in having a linear array of singly-coordinated sites which accept hydrogen bonds from triply-coordinated 
surface hydroxyls. The main differences are that the 0-0 distances are closer, the linear arrays are more closely spaced 
in the Y direction (0.518 nm vs 1.139 nm) and the hydrogen bonds donated by the triply coordinated hydroxyls are 
weaker. The hematite (012) surface should therefore prefer oxyanions with smaller 0-0 distances and have a higher 
overall specific adsorptive capacity than goethite (110). Assume for the sake of argument that binding is bidentate 
(eg. (=FeOH) 2 (P04) groups). For hematite (012) there are 2 bidentate sites per unit cell (0.288 nm 2 ) or 6.9 bidentate 
sites/nm 2 , whereas for goethite there is one bidentate site per unit cell (.358 nm 2 ), or 2.8 bidentate sites/nm 2 . Also, 
because of the close spacing of the hydroxide rows, hematite might be less effective in binding elongated adsorbates. 



( 001 ) Hydroxylated Surface Structure 

For the (001) surface, the first step toward hydroxylation is the adsorption of a water molecule directly over 
the three-fold coordinated iron ion (Figure [?]). For the minimized configuration, this breaks the symmetry and a 
rectangular unit mesh results. The adsorbed water molecule sits 207 pm above the (formerly) three-fold coordinated 
iron ion. The coordinates of atoms in the unit cell of the slab for hydrated and hydroxylated surface are given in 



Table XVIII and XIX, respectively. Again, as shown in Table XVII hydration of the (001) surface has almost no effect 
on surface relaxation. 

Surface hydroxylation (Figure ||) takes place by dissociating the adsorbed water molecule and protonating 1/3 of 
the three-fold coordinated surface oxygens. As for the (012) surface, hydroxylation has a significant effect on the 
surface relaxation. The top iron ion moves upward relative to the first oxygen layer increasing the distance between 



the top Fe and O planes to approximately the bulk value value. The hydroxide layer is not indicated in Table XVII 
moves about 22 pm above the oxygen ions not receiving protons in the first structural oxygen layer. There are no 
significant hydrogen bonding interactions at the surface. 

The surface hydroxyl speciation for (001) is two FesOH and two FeOH sites per unit cell (0.465 nm 2 ), resulting in 
a density of 4.3 sites/nm 2 of each type (seperately). There is no possibility for bidentate binding at the (001) surface, 
at least for the model arrived at here. 



Hydroxylation Energetics 

The energy of adsorption of water to the (012) surface in the configuration shown in Figure [| is -1.6 J/m 2 (the 
adsorbed water molecule is assumed to be removed from bulk ice at K). If, in addition, the water is allowed to 
dissociatively adsorb (as shown in Figure ^), the adsorption energy is -2.7 J/m 2 . Therefore, one obtains more energy 
for adsorption of the water layer than the 2 J/m 2 required to fracture the crystal. This is clearly not correct, since 
hematite crystals are stable in the presence of water. 

The same considerations apply to the "gas phase" hydroxylation of the (001) surface. Here the energy of adsorption 
of undissociated water in the configuration shown in Figure is -1.2 J/m 2 . The energy of dissociative adsorption in 
the configuration shown in Figure |] is -1.8 J/m 2 . Again the energy of dissociative adsorption exceeds the surface 
fracture energy (1.6 J/m 2 ). 

The surface enthalpy of hematite has been estimated in reference 0] by measuring the heat of dissolution of 
hematite crystals as a function of surface area p9| . A value of +0.770 J/m 2 was estimated for the surface enthalpy at 
70°C. Taking this result at face value, and assuming the experimental value is representative of the (012) surface, there 
is a discrepency of 1.4 J/m 2 between the experimental value and the theoretical estimate. The failure to obtain the 
correct sign could be due to either overestimating the binding energy of water, hydroxide, and protons to the surface, 
or underestimating the energy required to fracture the surface. As discussed above, nearly all theoretical estimates of 
the fracture energy are in fair agreement predicting that the surface energy is 1.6-2.0 J/m 2 yet experimental estimates 
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range from 6-40 J/m 2 . The discrepency of 1.4 J/m 2 is small in relation to this range of experimental values for the 
fracture energy. The same arguments apply if the predominant surface in the surface enthalpy experiments is (001). 
Meaningful discussion of the hydroxylation issue must wait until the overall discrepency between the theoretical and 
experimental fracture surface energies is better understood. We will, however, point out that our model predicts a 
greater amount of dissociation on (012) than (001), both in terms of "gas phase" calculations and qualitatively in the 
MD simulations of the hematite-water interface. 



VII. SUMMARY AND CONCLUSIONS 

We have presented a simple yet rigorous derivation of the 2-D Ewald summation method of Rhee et al. J8| and its 
generalization to polarizable systems. This provides a viable means of investigating the surfaces of oxide minerals 
using thousands of atoms. As one can see from the presented results of benchmarks against exact methods, the 
method of Rhee and co-workers is adequate for systems with aspect ratios {{z ma x — z m in)/min(L x ,L y )) less than 
approximately 0.5. We increased the order of the Taylor expansion for V m from 4 (as in |8|] and |9|]) and found this 
to yield significantly better accuracy. These expressions are given in the Appendix. 

We have applied the method to investigate the structures of hydroxylated and non-hydroxylated (012) and (001) 
surfaces of hematite (a-F^Os) in slab geometry, and solidly document the convergence of our results with respect 
to slab thickness. Our nonhydroxylated structures and energies are in good agreement with previous calculations on 
isostructural ev-A^Os. We used the flexibility and economy of our parameterized analytical potential approach to 



run simulations of solvated interfaces. These simulations revealed nonuniform hydroxylation as suggested in 35 1, but 
revealed end-member structures which were used in more precise "gas-phase" calculations. These calculations revealed 
that the effect of hydration on surface structure is small, but that hydroxylation has a significant effect. Hydroxylation 
appears to be energetically favored both in the gas-phase calculations and the solution phase simulations, with the 
(012) surface exhibiting a greater degree of dissociation of adsorbed water. The (001) surface has a complex solvated 
structure. Our simulations reveal the presence of both four-fold coordinated and six-fold coordinated iron ions at the 
solvated surface, with four-fold coordination predominant. Conversion of four-fold to six-fold coordination requires 
significant movement of top-layer iron atoms out into the aqueous solvent. This movement should be detected in x-ray 
reflectivity experiments designed to measure the location of the top layer of iron ions and their bound hydroxyl/aquo 
ligands in the direction normal to the surface plane [ j50|]5~if| . It is hoped that the present study will encourage 
experimental studies of hydroxylated hematite surfaces. 

Quantum mechanical simulations are potentially superior to those with parameterized analytical potentials because 
electronic structure is accounted for explicitly. However, these methods must always compromise some accuracy for the 
sake of feasibility. This means that even demonstrating convergence with respect to basis set size, densit y of fc -points 



mesh for Brillouin zone integration, system size requires tremendous effort. As shown in Tables XI and Kill, results 
of electronic structure based methods are not always in agreement. We therefore wish to emphasize the usefulness of 
parameterized analytical potential methods in extending the scope of theoretical calculations on mineral surfaces. 
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APPENDIX: MULTIPOLE EXPANSIONS 



The first, second, and third derivatives of V m in equatio n |2q are required for computation of energies and forces 
for charge-dipole and dipole-dipole interactions in equations |l| and |[ 
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FIG. 1. Nonhydrated (012) surface, (a) slab cross section of the system consisting of 3 layers of neutral O-Fe-O-Fe-O 
sequences. Large red atoms are oxygen, small orange atoms are iron, (b-d) closeup views of top O-Fe-O-Fe-O layer plus the 
top O layer from the middle neutral stacking sequence. Atoms are color coded according to depth. The top O layer is red, 
the middle O layer is magenta, and the bottom O layer is pink. The top O layer of the middle stacking sequence is grey. The 
top five-fold coordinated iron ions are brown, and the lower six-fold coordinated iron ions are orange. In the map view (d) the 
zig-zag oxygen ridges are emphasized with black vectors. 



FIG. 2. Nonhydrated (001) surface, (a) slab cross section of the system consisting of four layers of neutral Fe-O-Fe sequences. 
Large red atoms are oxygen, small orange atoms are iron, (b-d) closeup views of top two Fe-O-Fe layers but minus the bottom 
Fe layer from the 2nd neutral Fe-O-Fe stacking sequence. Atoms are color coded according to depth. The top O layer is red, 
and the bottom O layer is deep pink. The top three-fold coordinated ion ions are brown, the second iron layer is orange and 
the third iron layer is grey. 



FIG. 3. A snapshot (after 50 picoseconds) of the solvated (012) surface. The initial conditions of the simulation consisted 
were obtained by placing a 1.5 nm layer of water from a bulk water simulation over the surface shown in Figure [j](a). Atom 
types are the same as in Figure [l] except that both surface water molecules and solvent hydroxide ions are indicated by large 
yellow oxygen atoms. A water molecule is considered a surface water molecule if it is bound to at least 1 iron atom, and a 
hydroxide ion is considered solvated if it is not bound to any iron atoms. 



FIG. 4. A snapshot (after 50 picoseconds) of the solvated (001) surface. The initial conditions of the simulation consisted 
were obtained by placing a 1.5 nm layer of water from a bulk water simulation over the surface shown in Figure ^(a). Atom 
types are the same as in Figure ^| except that both surface water molecules and solvent hydroxide ions are indicated by large 
yellow oxygen atoms. A water molecule is considered a surface water molecule if it is bound to at least 1 iron atom, and a 
hydroxide ion is considered solvated if it is not bound to any iron atoms. 



FIG. 5. Three orthogonal views of the hydrated (012) surface. Atom types are the same as in Figure |l| except that adsorbed 
water molecules are indicated by the yellow oxygen atoms and protons are the small white atoms. Idealized periodic surface 
structures are those suggested from the MD simulations shown in Figure 



FIG. 6. Three orthogonal views of the hydroxylated (012) surface. Atom types are the same as in Figure [l] except that 
adsorbed hydroxide ions are indicated by the yellow oxygen atoms and protons are the small white atoms. Idealized periodic 
surface structures are those suggested from the MD simulations shown in Figure 



FIG. 7. Three orthogonal views of the hydrated (001) surface. Atom types are the same as in Figure g except that adsorbed 
water molecules are indicated by the yellow oxygen atoms and protons are the small white atoms. Idealized periodic surface 
structures are those suggested from the MD simulations shown in Figure 



FIG. 8. Three orthogonal views of the hydroxylated (001) surface. Atom types are the same as in Figure |l| except that 
adsorbed hydroxide ions are indicated by the yellow oxygen atoms and protons are the small white atoms. Idealized periodic 
surface structures are those suggested from the MD simulations shown in Figure H 
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TABLE I. The energy (in units of e 2 /L) and forces (in units of e 2 /L 2 ), between a charge +e at (0, 0, 0) and — e at (0.1, 0.1, 
0.1) L in an infinitely periodic 2-D array of square unit cells 



Method U F x 
Leckner (l max = 10) -5.77211 

Farry (k = oj -o.i iZLZ 18.8136 
This study -5.77184 18.8136 


Fy 

18.8136 
18.8136 


F z 

OA 1 COO 

zU. 1625 

20.1627 


TABLE II. The energy (in units of e 2 /L) and forces (in units of e 2 /L 2 ), between a 
-0.15, -0.2)L in an infinitely periodic 2-D array of square unit cells 


charge +e at (0, 0, 0) and 


-e at (-0.25, 


Method U F x 
Leckner (l max = 10) -2.82156 -4.60495 
Farry (n = o) -Z.olloi -4.6U495 
This study (7 = 0.1L short ) -2.8214336 -4.60323 


Fy 

-2.91052 

-2.911)02 

-2.909703 


F z 
-6.44748 
-6.44/4/0 

-6.446037 


TABLE III. The energy (in units of e 2 / ' L) and forces (in units of e 2 / L 2 ), between a 
-0.15, -0.5)L in an infinitely periodic 2-D array of square unit cells 


charge +e at (0, 0, 0) and 


-e at (-0.25, 


Method U F x 
Parry's -0.8061173 -0.6608550 
This study -0.8067679 -0.650034 
Error, % 0.08 1.64 


Fy 

-0.4514618 
-0.4461202 
1.18 


F z 

-6.723419 
-6.5550016 
2.50 



TABLE IV. The energy (in units of e 2 /A) and forces (in units of e 2 / A 2 ), between a charge +e at (0, 0, 0) and — e at (-0.25, 
-0.15, -0.1)L and also dipole fi = (—0.33, —0.16, 0.77) at the position of — e in an infinitely periodic 2-D array of square unit 
cells with the side length L = 10 A 



Method 


U 


F x x 10 2 


F y x 10 2 


F z x 10 2 


Parry's 


-0.33385337 


-7.280129828 


-4.56935990 


-8.077974417 


This study 


-0.336207073 


-7.2792028 


-4.56935929 


-8.077045 



TABLE V. The energy (in units of e 2 /A) and force (in units of e 2 / A 2 ) in a system consisting of a charge +2e at (0,0,0), 
a charge — e and dipole p = (—0.33,-0.16,0.77) at (—0.25,-0.15,-0.1)1/ and a charge -e and dipole fl = (0.3,0.6,0.7) at 
(0.5, 0.33, —0.05)1/ in an infinitely periodic 2-D array of square unit cells with the side length L — 10 A 



Method 


U 


F x 


Fy 


F z a 


Parry's 


-0.8925432834 


-0.13797318 


-0.062559269 


-0.2051223 


This study 


-0.89619000 


-0.13788360 


-0.0625851709 


-0.20480698 



a Only the force acting on charge +2e is given 
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TABLE VI. Slab unit cell coordinates for the ideal vacuum-terminated (012) surface. X, Y, and Z coordinates in A. 
Orthogonal unit mesh translation vectors are 5.564 A in X and 5.176 A in Y. 



atom 


X 


y 


z 


Ol 


0.374 


0.254 


1.169 


02 


2.689 


2.746 


4.819 


03 


3.605 


2.856 


-2.573 


04 


2.782 


0.000 


-3.870 


05 


2.066 


0.000 


3.868 


06 


5.206 


2.608 


-0.001 


07 


2.160 


2.746 


-4.820 


08 


1.244 


2.857 


2.572 


09 


4.474 


0.254 


-1.170 


O10 


1.692 


4.432 


-1.171 


Oil 


4.942 


1.940 


-4.820 


012 


4.026 


1.830 


2.572 


013 


4.848 


4.687 


3.868 


014 


2.424 


2.078 


-0.001 


015 


0.000 


4.686 


-3.870 


016 


0.823 


1.830 


-2.573 


017 


5.471 


1.941 


4.819 


018 


3.156 


4.432 


1.169 


Fel 


0.901 


1.056 


-0.786 


Fe2 


3.962 


3.526 


-4.531 


Fe3 


3.332 


3.656 


3.028 


Fe4 


1.516 


3.655 


-3.029 


Fe5 


0.887 


3.526 


4.530 


Fe6 


3.947 


1.056 


0.785 


Fe7 


3.669 


1.161 


4.530 


Fe8 


1.165 


3.630 


0.785 


Fe9 


4.298 


1.030 


-3.029 


FelO 


1.180 


1.160 


-4.531 


Fell 


0.550 


1.031 


3.028 


Fel2 


3.683 


3.630 


-0.786 



TABLE VII. Changes in interplanar spacing due to relaxation of hematite (012) surface as a function of size of the simulation 
cell 

Change from bulk value 
Number of Fe2 03 molecules 



value (pm) 


144 


216 


288 


35 


-24% 


-18% 


-18% 


80 


-24% 


-17% 


-18% 


80 


8% 


5% 


5% 


35 


49% 


32% 


34% 


149 


-15% 


-6% 


-7% 



1G 



TABLE VIII. Surface energy 7 (J/m 2 ) for non-hydroxylated (012) face of hematite calculated using the 3-D and 2-D periodic 
boundary conditions 



Number of Fe2C>3 molecules 73-D relaxed 72-D relaxed 73-D unrelaxed 

72 2.335 2.325 2.903 

144 1.882 1.874 2.757 

216 1.994 2.022 2.775 

288 1.999 2.094 2.787 



TABLE IX. Slab unit cell coordinates for the ideal vacuum-terminated (001) surface. X, Y, and Z coordinates in A. 
Orthogonal unit mesh translation vectors are 8.976 A in X and 5.176 A in Y. 



atom 


X 


y 


z 


Ol 


0.000 


4.967 


-1.188 


02 


1.410 


2.409 


3.428 


03 


4.488 


2.380 


-1.188 


04 


5.898 


4.997 


3.428 


05 


7.360 


2.497 


3.428 


06 


5.847 


0.000 


-1.188 


07 


1.359 


2.587 


-1.188 


08 


2.872 


5.084 


3.428 


09 


8.798 


5.011 


3.428 


O10 


7.231 


2.365 


-1.188 


Oil 


2.743 


4.952 


-1.188 


012 


4.310 


2.423 


3.428 


013 


4.317 


0.752 


-3.428 


014 


5.727 


3.369 


1.188 


015 


8.805 


3.340 


-3.428 


016 


1.239 


0.782 


1.188 


017 


8.856 


3.162 


1.188 


018 


7.343 


0.665 


-3.428 


019 


2.855 


3.252 


-3.427 


O20 


4.368 


0.574 


1.188 


021 


5.905 


3.326 


-3.428 


022 


7.472 


0.796 


1.188 


023 


1.418 


0.738 


-3.428 


024 


2.984 


3.384 


1.188 


Fel 


4.359 


4.167 


-3.863 


Fe2 


5.856 


1.580 


0.290 


Fe3 


8.847 


1.580 


-3.863 


Fe4 


1.368 


4.168 


0.290 


Fe5 


5.855 


1.580 


-2.603 


Fe6 


7.352 


4.168 


2.209 


Fe7 


1.367 


4.168 


-2.602 


Fe8 


2.864 


1.581 


2.209 


Fe9 


8.847 


1.581 


-0.289 


FelO 


5.856 


1.582 


3.864 


Fell 


1.368 


4.169 


3.864 


Fel2 


4.359 


4.169 


-0.289 


Fel3 


8.848 


1.581 


2.603 


Fel4 


7.351 


4.168 


-2.209 


Fel5 


2.864 


1.581 


-2.209 


Fel6 


4.360 


4.168 


2.603 
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TABLE X. Changes in interplanar spacing due to relaxation of hematite (001) surface as a function of size of the simulation 
cell 



Change from bulk value 
Number of Fe2 03 molecules 



Bulk value (pm) 


168 


224 


280 


336 


86 


-53% 


-50% 


-49% 


-49% 


86 


-3% 


-3% 


-2% 


-3% 


63 


-43% 


-38% 


-43% 


-41% 


86 


24% 


19% 


21% 


21% 



TABLE XL Surface Energies of Non-hydroxylated Fe/Al Sesquioxide (J/m 2 ): Comparison with previous electronic structure 
calculations 

a-Fe2 03 a-Ai2 03 

Ours LDA H HF @ HF ||] 

(001) 

Unrelaxed 4.24 3.77 3.3 5.0 

Relaxed 1.64 1.76 2.0 

(012) 

Unrelaxed 2.79 2.95 

Relaxed 2.00 



TABLE XII. Surface energy 7 (J/m 2 ) for non-hydroxylated (001) face of hematite calculated using the 3-D and 2-D periodic 
boundary conditions 



Number of Fe2 03 molecules 
112 
168 
224 
280 
336 



73 



-n relaxed 
1.63 
1.54 
1.64 
1.64 
1.65 



72 -d relaxed 
1.62 
1.46 
1.63 
1.64 
1.60 



73 



-d unrelaxed 
4.32 
4.22 
4.22 
4.23 
4.24 



TABLE XIII. Nonhydroxylated (001) surface. Percent relaxation in z direction 





Ours 


LDA y 


TBTE Q 


hf y 


hf y 


Fe/Al 










-49 


-86 


-87 


-50 


-50 

















-3 


+3 


-33 


-5 




Fe/Al 


-41 


-54 


+20 


-8 




Fe/Al 


21 


25 









O 
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TABLE XIV. Slab unit cell coordinates for the hydrated (012) surface. X, Y, and Z coordinates in A. Orthogonal unit mesh 
translation vectors are 5.564 A in X and 5.176 A in Y. 



atom 


X 


y 


z 


HI 


1.856 


0.028 


-6.877 


H2 


1.122 


1.264 


-7.312 


H3 


3.728 


1.264 


7.314 


H4 


2.987 


0.033 


6.879 


H5 


0.946 


3.942 


7.312 


H6 


0.221 


0.000 


6.883 


H7 


4.628 


0.003 


-6.882 


H8 


3.906 


3.938 


-7.313 


Ol 


1.135 


0.643 


-6.531 


02 


0.377 


0.519 


1.166 


03 


2.721 


3.024 


4.809 


04 


3.596 


3.123 


-2.575 


05 


2.782 


0.264 


-3.846 


06 


2.068 


0.264 


3.846 


07 


5.207 


2.873 


0.001 


08 


3.712 


0.644 


6.533 


09 


2.129 


3.024 


-4.809 


O10 


1.253 


3.123 


2.576 


Oil 


4.472 


0.520 


-1.165 


012 


0.931 


4.563 


6.532 


013 


1.691 


4.690 


-1.165 


014 


4.910 


2.185 


-4.809 


015 


4.035 


2.086 


2.576 


016 


4.849 


4.945 


3.847 


017 


2.425 


2.337 


0.000 


018 


0.000 


4.945 


-3.845 


019 


0.815 


2.086 


-2.575 


O20 


5.503 


2.186 


4.810 


021 


3.923 


4.560 


-6.533 


022 


3.159 


4.690 


1.166 


Fel 


0.902 


1.315 


-0.786 


Fe2 


3.957 


3.813 


-4.546 


Fe3 


3.332 


3.908 


3.022 


Fe4 


1.518 


3.908 


-3.021 


Fe5 


0.893 


3.813 


4.547 


Fe6 


3.948 


1.315 


0.787 


Fe7 


3.675 


1.396 


4.548 


Fe8 


1.166 


3.894 


0.787 


Fe9 


4.300 


1.302 


-3.021 


FelO 


1.175 


1.396 


-4.546 


Fell 


0.550 


1.301 


3.021 


Fel2 


3.684 


3.894 


-0.786 
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TABLE XV. Slab unit cell coordinates for the hydroxylated (012) surface. X, Y, and Z coordinates in A. Orthogonal unit 
mesh translation vectors are 5.564 A in X and 5.176 A in Y. 



atom 


X 


y 


z 


HI 


1.763 


0.441 


-7.029 


H2 


2.936 


3.749 


5.723 


H3 


2.782 


0.442 


7.030 


H4 


1.609 


3.748 


-5.723 


H5 


0.000 


0.068 


7.030 


H6 


4.391 


1.937 


-5.723 


H7 


0.154 


1.938 


5.723 


H8 


4.545 


0.069 


-7.029 


Ol 


1.070 


0.429 


-6.342 


02 


0.218 


0.764 


1.144 


03 


2.587 


3.251 


4.895 


04 


3.364 


3.334 


-2.612 


05 


2.663 


0.510 


-3.766 


06 


1.881 


0.510 


3.767 


07 


5.054 


3.121 


0.000 


08 


3.475 


0.430 


6.343 


09 


1.958 


3.251 


-4.895 


O10 


1.181 


3.334 


2.612 


Oil 


4.327 


0.764 


-1.143 


012 


0.693 


0.080 


6.343 


013 


1.545 


4.922 


-1.143 


014 


4.740 


2.434 


-4.895 


015 


3.963 


2.352 


2.612 


016 


4.663 


5.176 


3.767 


017 


2.272 


2.564 


0.000 


018 


5.445 


0.000 


-3.767 


019 


0.582 


2.352 


-2.612 


O20 


5.369 


2.434 


4.895 


021 


3.852 


0.080 


-6.342 


022 


3.000 


4.922 


1.143 


Fel 


0.749 


1.546 


-0.802 


Fe2 


3.912 


4.265 


-4.733 


Fe3 


3.149 


4.113 


2.971 


Fe4 


1.396 


4.113 


-2.970 


Fe5 


0.633 


4.265 


4.733 


Fe6 


3.796 


1.546 


0.802 


Fe7 


3.415 


1.421 


4.733 


Fe8 


1.014 


4.140 


0.802 


Fe9 


4.178 


1.573 


-2.970 


FelO 


1.130 


1.421 


-4.733 


Fell 


0.367 


1.573 


2.971 


Fel2 


3.531 


4.140 


-0.802 
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TABLE XVI. Interplanar spacings (picometers) of hydrated and hydroxylated hematite (012) (as shown in Figure ^|) com- 
pared with bulk and relaxed nonhydroxylated surface structures 



O 
Fe 

o 

Fe 
O 
O 



bulk 

35 

80 

80 

35 
149 



relaxed 
29 (-18%) 
66 (-17%) 

84 (5%) 
46 (32%) 
140 (-6%) 



hydrated 
26 (-26%) 
70 (-12%) 
82 (3%) 
45 (29%) 
141 (-5%) 



hydroxylated 
16 (-54 %) 
97 (21 %) 
80 (-1 %) 
36 (4 %) 
147 (-1 %) 



TABLE XVII. Interplanar spacings (picometers) of hydrated and hydroxylated hematite (001) (as shown in Figure (?]) 
compared with bulk and relaxed nonhydroxylated surface structures 



Fe 

o 

Fe 
Fe 

o 



bulk 

86 
86 
63 
86 



relaxed 
44 (-49%) 
83 (-4%) 
39 (-37%) 
102 (19%) 



hydrated 
48 (-44%) 

83 (-3%) 
38 (-39%) 
102 (19%) 



hydroxylated 
89 (4%) 
78 (-9%) 
42 (-32%) 
98 (14%) 
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TABLE XVIII. Slab unit cell coordinates for the hydrated (001) surface. X, Y, and Z coordinates in A. Orthogonal unit 
mesh translation vectors are 8.976 A in X and 5.176 A in Y. 



atom 


X 


y 


z 


HI 


4.989000 


3.646000 


-6.607000 


H2 


3.927000 


4.686000 


-6.662000 


H3 


0.400000 


2.017000 


-6.664000 


H4 


8.314000 


0.966000 


-6.607000 


H5 


6.296000 


1.119000 


6.663000 


H6 


5.210000 


2.132000 


6.584000 


H7 


1.875000 


4.760000 


6.608000 


H8 


0.795000 


3.739000 


6.640000 


Ol 


0.000000 


4.964000 


-1.193000 


02 


1.410000 


2.409000 


3.422000 


03 


4.490000 


2.376000 


-1.190000 


04 


5.904000 


4.997000 


3.419000 


05 


7.366000 


2.494000 


3.416000 


06 


5.851000 


0.000000 


-1.194000 


07 


1.364000 


2.590000 


-1.189000 


08 


2.877000 


5.086000 


3.426000 


09 


8.800000 


5.009000 


3.407000 


O10 


7.231000 


2.367000 


-1.197000 


Oil 


2.743000 


4.954000 


-1.186000 


012 


4.314000 


2.427000 


3.423000 


013 


4.316000 


0.753000 


-3.429000 


014 


5.731000 


3.373000 


1.179000 


015 


8.810000 


3.340000 


-3.434000 


016 


1.245000 


0.785000 


1.181000 


017 


8.857000 


3.159000 


1.177000 


018 


7.344000 


0.663000 


-3.435000 


019 


2.855000 


3.255000 


-3.428000 


O20 


4.369000 


0.574000 


1.183000 


021 


5.906000 


3.322000 


-3.435000 


022 


7.477000 


0.795000 


1.174000 


023 


1.420000 


0.740000 


-3.421000 


024 


2.990000 


3.383000 


1.185000 


025 


4.447000 


4.188000 


-5.980000 


026 


8.867000 


1.499000 


-5.979000 


027 


5.779000 


1.604000 


5.969000 


028 


1.324000 


4.241000 


5.969000 


Fel 


4.359000 


4.168000 


-3.910000 


Fe2 


5.859000 


1.581000 


0.293000 


Fe3 


8.851000 


1.579000 


-3.909000 


Fe4 


1.372000 


4.168000 


0.295000 


Fe5 


5.856000 


1.579000 


-2.601000 


Fe6 


7.356000 


4.167000 


2.199000 


Fe7 


1.368000 


4.169000 


-2.596000 


Fe8 


2.867000 


1.582000 


2.205000 


Fe9 


8.848001 


1.581000 


-0.307000 


FelO 


5.861000 


1.582000 


3.898000 


Fell 


1.369000 


4.170000 


3.898000 


Fel2 


4.362000 


4.168000 


-0.304000 


Fel3 


8.853000 


1.581000 


2.584000 


Fel4 


7.353000 


4.167000 


-2.217000 


Fel5 


2.865000 


1.582000 


-2.211000 


Fel6 


4.365000 


4.171000 


2.590000 
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TABLE XIX. Slab unit cell coordinates for the hydroxylated (001) surface. X, Y, and Z coordinates in A. Orthogonal unit 
mesh translation vectors are 8.976 A in X and 5.176 A in Y. 



atom 


X 


y 


z 


HI 


8.478001 


5.121000 


4.469000 


H2 


3.991000 


2.534000 


4.470000 


H3 


6.468000 


3.444000 


-4.468000 


H4 


1.980000 


0.856000 


-4.469000 


H5 


4.417000 


4.322000 


-6.994000 


H6 


8.902000 


1.740000 


-6.994000 


H7 


6.046000 


1.651000 


6.994000 


H8 


1.547000 


4.245000 


6.994000 


Ol 


0.149000 


5.068000 


-1.204000 


02 


1.585000 


2.589000 


3.353000 


03 


4.638000 


2.480000 


-1.202000 


04 


6.073000 


0.000000 


3.355000 


05 


7.481000 


2.621000 


3.362000 


06 


5.992000 


0.073000 


-1.164000 


07 


1.503000 


2.660000 


-1.164000 


08 


2.993000 


0.033000 


3.363000 


09 


8.954000 


5.145000 


3.581000 


O10 


7.321000 


2.500000 


-1.156000 


Oil 


2.832000 


5.086000 


-1.154000 


012 


4.467000 


2.557000 


3.582000 


013 


4.386000 


0.800000 


-3.353000 


014 


5.819000 


3.495000 


1.203000 


015 


8.875000 


3.390000 


-3.355000 


016 


1.331000 


0.906000 


1.202000 


017 


8.955000 


3.317000 


1.163000 


018 


7.465000 


0.769000 


-3.362000 


019 


2.977000 


3.356000 


-3.362000 


O20 


4.467000 


0.729000 


1.165000 


021 


5.991000 


3.418000 


-3.581000 


022 


7.626000 


0.891000 


1.155000 


023 


1.503000 


0.830000 


-3.582000 


024 


3.137000 


3.476000 


1.157000 


025 


4.488000 


4.316000 


-6.081000 


026 


0.000000 


1.726000 


-6.081000 


027 


5.974000 


1.661000 


6.081000 


028 


1.482000 


4.249000 


6.080000 


Fel 


4.317000 


4.316000 


-4.252000 


Fe2 


5.988000 


1.706000 


0.358000 


Fe3 


8.805000 


1.729000 


-4.252000 


Fe4 


1.500000 


4.292000 


0.357000 


Fe5 


5.945000 


1.611000 


-2.577000 


Fe6 


7.426000 


4.272000 


2.153000 


Fe7 


1.457000 


4.197000 


-2.577000 


Fe8 


2.938000 


1.683000 


2.155000 


Fe9 


8.959000 


1.685000 


-0.359000 


FelO 


6.142000 


1.661000 


4.252000 


Fell 


1.655000 


4.249000 


4.250000 


Fel2 


4.471000 


4.272000 


-0.357000 


Fel3 


0.024000 


1.779000 


2.579000 


Fel4 


7.520000 


4.294000 


-2.152000 


Fel5 


3.030000 


1.703000 


-2.152000 


Fel6 


4.512000 


4.365000 


2.579000 
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